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ABSTRACT 



We examine the nonlinear magnetofrictional extrapolation scheme using the solar active region model by Titov and Demoulin as test 
field. This model consists of an arched, line-tied current channel held in force-free equilibrium by the potential field of a bipolar flux 
distribution in the bottom boundary. A modified version, having a parabolic current density profile, is employed here. We find that the 
equilibrium is reconstructed with very high accuracy in a representative range of parameter space, using only the vector field in the 
bottom boundary as input. Structural features formed in the interface between the flux rope and the surrounding arcade — "hyperbolic 
flux tube" and "bald patch separatrix surface" — are reliably reproduced, as are the flux rope twist and the energy and helicity of the 
configuration. This demonstrates that force-free fields containing these basic structural elements of solar active regions can be obtained 
by extrapolation. The influence of the chosen initial condition on the accuracy of reconstruction is also addressed, confirming that the 
initial field that best matches the external potential field of the model quite naturally leads to the best reconstruction. Extrapolating 
the magnetogram of a Titov-Demoulin equilibrium in the unstable range of parameter space yields a sequence of two opposing 
evolutionary phases which clearly indicate the unstable nature of the configuration: a partial buildup of the flux rope with rising free 
energy is followed by destruction of the rope, losing most of the free energy. 

Key words. Magnetic fields — Sun: corona — Sun: Surface magnetism 



1. Introduction 

The computation of the coronal magnetic field from bound- 
ary data is the only available technique to obtain fully three- 
dimensional information, i.e., the complete field vector in the 
whole coronal volume bounded below by the magnetogram. 
Deductions of the field from spectral lines formed in the coronal 
plasma, from radio maps across a range of microwave frequen- 
cies, or from observations of specific structures such as promi- 
nences, are restricted to two-dimensional projections or to cuts 
through the volume, or do not yield the full vector. The three- 
dimensional information is required to advance a wide array of 
issues in coronal physics, for example the structure of promi- 
nences, the onset of eruptions, and how the coronal field is con- 
nected to the Sun's interior and the solar wind. 

In equilibrium, the coronal magnetic field is nearly force-free 
in the majo r part of active regions, except their upper peripher y 
(see, e.g., iMetcalf et al.lll995l: iMoon et al.ll2QQ2t [G^l2:00lh . 
Hence, to a very good approximation it can be described by 



VxB = aB 



with 



V-B = 0. 



(1) 



Here a(x) is a scalar function that is constant along each indi- 
vidual field line (which follows directly from Eqs. |[T|) but in 
general has diflTerent values on diflTerent field lines. 

The computation of force-free coronal fields can be formu- 
lated as a boundary- value problem for Eqs. ([T]). Such extrap- 
olation can be performed at diflTerent levels of approximation 
for the field. If a is assumed to be constant, then the first of 
Eqs. ([T]) is linear, and only the knowledge of the normal field 



component is needed as boundary condition for the extrapola- 
ti on. The solution of the linear problem was given in c l osed form 
in 'Nakagawa & Raadu' ('I972h: IChiu & HiltonI ([T977h : ISeehafed 
(1978). Linear methods have several intrinsic limitations, the 
most severe one being that they cannot match most observed 
magnetograms closely becau se a is often found to vary strongly 
across active regions (e.g., iRegnier et al.l 120021) . The particu- 
lar case of a = yields a potential field, which is inadequate 
for many purposes, since it does not contain any free energy 
to power coronal activity. Therefore, it is necessary to proceed 
beyond the linear approximation and to solve the nonlinear ex- 
trapolation problem, permitting Qr(jt:) to vary across the field 
(ISakurailll98il:lMcClvmont et al.lll997l:rAmari et al.lll999h . Due 
to the nonlinearity, such extrapolation requires a numerical ap- 
proach. Various numerical schemes and implementations have 
been deve loped to construct a solution; for a recent review see 
ISchriiver et al. (2006). 

A pa rtially alternative stra t egy, the flux rope insertio n 
method (van Ballegooijen 2004; 'van Ballegooije n et al.l l2007h . 
permits to model structures that are supposed to contain a flux 
rope, particularly filament channels. Here an extrapolated poten- 
tial field is replaced by a flux rope in part of the volume, and the 
resulting configuration is numerically relaxed to find an equilib- 
rium. The parameters of the inserted rope are varied until the 
field lines of the resulting equilibrium match observed features, 
e.g., threads in a filament or overlying loops. 

A convenient tool for both nonlinear force-free extrapolation 
and the flux rope insertion method is the so-called magnetofric- 
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tio nal relaxation of the field to force-free numerical equilib- 
ria ("Chodura & Schluter"1981'; 'Yang et alJll98S ICraig & Snev3 
1986 : Roumeliotis 1996). In a previous paper dValori et al. 
20051) . we described our implementation of magnetofrictional 
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nonlinear extrapolation and its application to reconstruct a rel- 
atively simple force-free magnetic field, which contains an ap- 
proximately current-neutralized flux rope rooted in the flux con- 
centrations of a simple bipolar magnetogram (Torok & Kliem 
l2003l) . That paper included a thorough study of the dependence 
of the reco n structi on quality on the numerical resolution. In 
IValori et aP (l2007l) an improved implementation of the mag- 
netofrictional method allowed us to significantly increase the ac- 
curacy i n the reconstruction of the well known force-free equi- 



libria by lLow & Lou (Il990|) in comparison to similar extrapola- 
tions summarized in lSchrijver et aiT(l2006h . Applications of this 



code to measured magnetogra ms are included in Metcal f et al.l 
(120081) . ISchriiver et al. I (12008 ah . and lDeRosa etai.1 (I2009h . 

Here we extend the testing of our magnetofrictional code 
in the reconstruction of solar-relevant solutions of Eqs. ([B 
by applying it to the fl ux rope equilibrium constructed by 
Iritov & DemoulinI (Il999l hereafter TD). This configuration de- 
scribes structures that emerge fro m the solar interior already 
twisted and carry a net current (e.g.. lLites et al.lll995l: lLeka et al.l 
[1996: Wheatland 20 00), and it is in this sense complemen- 
tary to the one in Valori et al.l (|2005[) which represents the 
case of an originally potential field twisted solely by photo- 
spheric displacements. The relevance of the TD equilibrium 
as a model for solar active regions and eruptive configurations 
was demonstrated by many inves tigations (e.g.. Zarro et al. 
20001: iKliemetaLl 120 04: Torok & Kfieml |2Q05|; ISchriiver et aL 
2008bl: iMcKenzie & Canfieldii2008h . In particular, it has been 



shown t hat the TD flux rope can be subject to ideal MHD insta- 
bilities (^Roussev et al.' 2003; Torok et al. 20041; iTorok & Kliem 
120071: llsenberg & Fo rbes 2007). Depending on its parameters, 
the equilibrium can include a hyperbolic flux tube (HFT) or bald 
patches (BPs) (see Sect. [2l for detail). Such structural features 
are thought to play an important role in the initiation and evo- 
lution of coronal mass ejections (CMEs) and their associated 
flares. Moreover, the occurrence of BPs changes the topology 
of the field. 

Our goal in this paper is to provide a detailed analysis of 
the reconstruction capabilities of the magnetofrictional code for 
this test field. We address the flux rope morphology, including 
the twist and the structural features HFT and BP, the energy and 
helicity contents, and the influence of the initial condition used 
by the extrapolation code. We apply the code also to an unsta- 
ble case, expecting this to be useful in judging extrapolations 
of magnetograms that were taken during the initial phase of an 
eruption or immediately prior to its onset. Finally, we will com- 
pare our results to t he prev ious reconstruction of the TD field by 
IWiegelmann et al.l (l2006al) . which employed a diff'erent extrapo- 
lation scheme. 

Using a force-free test field implies that the magnetogram 
is compatible with the force-free hypothesis of the extrapola- 
tion method. Therefore, in the following we do not need to 
consider deviations from force-freeness in the photosphere, the 
ambiguity of the transverse field direction, noise, and insuf- 
ficient resolution, all of which are associated with measured 
magnetograms. Attemp ts to deal with some of th e se problems 
are di sc ussed, e.g., i n McClvmont et al. (1997), Amari et al. 
(Il997h . IWiegelmann et al. (2006 b), Fuhrmann et al. (2007) 
Wtcalf etal. (2008), Wiegelmann et al.l (l2008h . ISchriiver etaP 
(■2008a) . and .DeRosa et al.. ' (.2009.) . 
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Fig. 1. "Line-of- sight" magnetogram B^ix, y, 0) with superim- 
posed arrows representing the transverse field (top), and force- 
free parameter a(x, y, 0) (bottom), for the reference fields of the 
cases High_HFT (left) and BP (right). 



The paper is organized as follows. The TD equilibria used 
below are described in Sect. [2l the extrapolation method is 
briefly outlined in Sect. [3l and the results are presented and dis- 
cussed in Sect. IH Our conclusions are summarized in Sect. \5\ 



2. Test equilibria 

iTitov & DemoulinI (1 19991 TD) constructed a force-free equilib- 
rium of a toroidal current channel to obtain an analytical model 
of a solar active region that contains free magnetic energy. Being 
force free and carrying a net current, the equilibrium belongs to 
the cl ass of to kamak equilibria, also known as Shafranov equi- 
libria dShafrano v 1966), which require an external poloidal field, 
i.e., one due to sources other than the current channel, to balance 
the always outward directed Lorentz self-force (hoop force) of 
the bent current channel. TD used a pair of magnetic sources, 
placed symmetrically to the torus plane at the symmetry axis of 
the torus, to model the external poloidal field, ^ep, of a solar 
active region. They also added a line current running along the 
symmetry axis to include an external toroidal field, B^t, which 
permits to model the magnetic shear typically present in active 
regions and prevents the twist profile from becoming infinite at 
the surface of the torus. The torus is arranged vertically, with its 
symmetry axis submerged below the plane {z = 0}, which rep- 
resents the photosphere, resulting in a smooth, divergence-free 
field in the coronal volume {z > 0}. Inside the torus, the field is 
obtained by modifying the force-free field of a locally straight 
current channel of prescribed radial profile to match the external 
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potential field smoothly at the surface. The original version of 
the equilibrium uses a uniform radial profile of the current den- 
sity, resulting, after the matching, in a profile Jip) that increases 
moderately toward the surface at p = a, where it drops to zero 
[J{p > a) = 0] . The requirement of force-freeness implies that 
a toroidal field component is always present inside the torus, re- 
gardless of the value of 5et- The force-freeness within the current 
channel (i.e., the equilibrium condition) is satisfied only approx- 
imately, but numerical relaxation runs verify that the approxi- 
mation is very go od for high toroidal aspect ratio, Rja » 1 
(iKliem et al .1120041) , and tolerable down to mod erate and likely 
more r ealistic aspect ratios of order RIa ~ 2 ( Schrijver et al.l 
I2008bh . 

If characterized by its field lines, the equilibrium represents 
an arched, line-tied flux rope that is embedded in a (gener- 
ally sheared) potential field of arcade topology; the flux rope 
has a somewhat larger cross-section than the current channel. 
The "line-of- sight" magnetogram, B^ix, y, 0), possesses four flux 
concentrations, a pair of "sunspots" above the point sources, 
which are located at (±L, 0, -J), and a pair of "satellite polar- 
ities" in the intersection between the current channel and the 
"photosphere" (Fig. [T]). (Depending on the parameters chosen 
for the particular realization of the equilibrium, the satellite po- 
larities can actually encompass more flux than the sunspots.) The 
neutral line of the normal field component in the magnetogram, 
given by B^ix, y, 0) = 0, runs approximately parallel to the mag- 
netic flux rope axis beneath its apex and runs relatively straight 
between the two positive and the two negative flux concentra- 
tions, i.e., the "active region" is to be characterized as essentially 
bipolar (for a quadrupolar active region in the classical sense one 
would require the neutral line to meander substantially between 
opposite flux concentrations or even to be split, with one or sev- 
eral parts forming a closed loop). 

In a wide range of parameter space, approximately bounded 
by the conditions that d < R - a and ^et < ^ep under the 
rope, the interface between the flux rope and the surrounding 
arcade of field lines is a separatrix surface (or quasi- separatrix 
layer), a cross which the field line mapping varies abruptly (or 
rapidly) (iPriest & Demoulinlll995l:lDemoulin et al.| [T996). In the 
following we will collectively refer to these as separatrix. They 
introduce one or both of the following structural features, de- 
pending upon the parameters of the equilibrium. If the separa- 
trix intersects itself under the flux rope, a line of X-type topol- 
ogy results, which is a concentric ring of radius smaller than 
{R - a) lying in the torus plane [x = 0}. The flux bundle in the 
vicinity of t his line is often referred to as h yperbolic flux tube, 
HFT (.Titov eTaP I1999L |2002|: lTitovll2007h. The HFT pinches 
under external perturbations (iTitov et al.ll2003l : iGalsgaard et al.l 
'2003); thus it plays the role of a seed for the vertical flare 
current sheet in eruptions (Toroket al.ll2004l) and is important 
for the CME-flare relationship. If the separatrix touches the 
photosphere tangentially, it does so at the neutral line and at 
these locations the field points in the "inverse" direction (from 
the negative-polarity side of the neutral line to the positive- 
polarity side, see Fig. O. Such features are referred to as bald 
patch, BP ([Se ehafer 1986: lTitov et al.lll993l:lBungev et al.lll996l: 
iTitov & D^m oulin 1999i). They have been suggested to be sites 
where partial eruptions ori ginate and where sigmoidal coronal 
sources form (iGibson et al.ll2004l:lGreen & Kfiemll2009h . 

The bipolar structure of the external field implies the ab- 
sence of n ull points. These relevant topological features (e.g., 
lAntiochosi ri998) can be included by generalizing the external 
field, however, we leave this step for future work. 
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Fig. 2. Contours of B^{x,y,z - A), with superimposed arrows 
representing the transverse field in the vicinity of the neutral 
line of B^, for the reference (top) and reconstructed (bottom) BP 
cases. Note the inverse crossing of the neutral line, the charac- 
teristic signature of a BP, in the middle of the plots. 



With regard to the reconstruction task, it is important to note 
that TD equilibria which contain a separatrix are structurally 
quite diff'erent from any potential field that can be constructed 
from their corresponding line-of- sight magnetogram. The poten- 
tial field will not include a flux rope. In the majority of cases, 
HFT or BP(s) will not be present either, and in case they happen 
to exist, they will have very difl'erent locations than in the TD 
field. At the same time, these equilibria appear to be the sim- 
plest force-free configurations currently considered in coronal 
physics which are structurally (and in the presence of BP(s) even 
topologically) diff'erent from a sheared-arcade field. They con- 
sist of a single current channel and the minimum two flux sys- 
tems required: the flux in the rope and the external poloidal flux. 
These are the generic elements of more complex force-free con- 
figurations that carry a net current. The requirement of a struc- 
tural or even topological change makes the reconstruction quali- 
tatively diff'erent from reconstructions of Low & Lou fields or of 
fields obtained by shearing or twisting the photospheric flux with 
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smoo th profiles of substantial symmetry (e.g., iTorok & Klieml 
l2QQ3h . Therefore, the reconstruction of the TD field is a major 
and necessary step in verifying the capabilities of nonlinear ex- 
trapolation schemes. 

An impression of the difficulty faced by the extrapolation 
code can be obtained from the field line picture. Fig. [3] If an 
HFT is present, the flux connecting the satellite polarities must 
"penetrate" the flux connecting the sunspots in the course of 
the reconstruction. This is due to the fact that, in the poten- 
tial field that is used as initial condition for the extrapolation, 
the flux in the satellite polarities connects to the corresponding 
sunspot of opposite polarity, resulting in two main flux bundles 
lying side by side. In the final TD field, however, the connec- 
tion between the four flux concentrations in the magnetogram is 
like a cross (which can be traced back to the requirement that a 
non-neutralized current channel in force-free equilibrium must 
be stabilized by an external poloidal field). This results in mag- 
netic connections between the sunspots below and above the flux 
rope, i.e., the flux rope passes in between these two parts of the 
sunspot flux (this can be seen in Fig. [3] and, more clearly, in 
Fig. [8] below). If a BP is present, it is obvious that the field lines 
in and immediately above it must reverse the direction of neutral 
line crossing in proceeding from the initial potential field to the 
solution. This implies that field lines of the initial potential field 
passing low over the BP location, which are simple short loops 
rooted in the vicinity of the BP, must be replaced by field lines 
that are typically long and rooted at very diff'erent locations. 

We note that the configuration considered in iMetcalf et al.l 
(l20Q8h represents an intermediate case between purely sheared 
configurations (like, e.g., the family of Low & Lou fields) and 
the TD field. It was obtained by means of the flux rope insertion 
method such that the magnetic axis of the rope passes at a low 
height of only a few Mm abo ve the photosph ere along much of 
the length of the rope (van Ball egooijenl2Q04l) . This corresponds 
to a strongly submerged TD rope with R- a < d < R-\-a and con- 
tains a separatrix which m ay or may not form an HFT or a BP. 
The configuration used in iMetcalf et al.l (l2008h includes a BP, 
but the extrapolations presented in that paper were not checked 
for the reconstruction of this feature. 

We have searched published extrapolations of solar magne- 
togram data for the occurrence of BPs (inverse field direction 
at the neutral line) and of the flux penetration feature charac- 
teristic of an HFT There appear to be only two cases of a BP 
(ICanou et al .1120091: iGuo e t al. 2010) and one case of flux pene- 
tration (Reg nier & Am ari 2004). 

For the purposes of the present paper, the original TD equi- 
librium is modified in two ways. First, in place of the line cur- 
rent, we use a pa ir of dipoles as sources of the external toroidal 
field (iTorok & Kliem.200 5). This yields a compact distribution 
of magnetic flux in the magnetogram plane and a steeper de- 
crease of the external field with height than the original TD field, 
both of which are more realistic with regard to solar active re- 
gions. The dipoles are placed below the footpoints of the flux 
rope and point nearly vertically, such that the field lines con- 
necting them through the torus section in {z > 0} match the ge- 
ometry of the torus section closely. The conservation of toroidal 
flux in the rope now requires the minor radius of the torus to 
vary along its length. Since ^et still points approximately in the 
direction of the current, its strength remains a free parameter 
to a good approximation for parameters of interest here, similar 
to the original equilibrium. The break of the toroidal symmetry 
also causes the magnetic axis of the rope and the HFT to bulge 
out of the plane {x = 0}, but the efl'ect remains small, since ^et 
varies only moderately along the coronal flux rope for the cho- 



Table 1. Parameters of the TD equilibria 



Case 


R 


a 


d 


L 


^0 


m 


High_HFT 


1.83 


0.67 


0.83 


0.50 


-1.75 


-65.82 


Low_HFT 


1.83 


0.67 


0.83 


0.83 


-1.75 


-65.96 


No_HFT 


1.83 


0.67 


0.83 


1.17 


-1.75 


-65.80 


BP 


1.83 


0.90 


0.83 


1.17 


-1.75 


-38.51 


Unstable 


1.83 


0.41 


0.83 


0.83 


5.00 


-186.7 



Normalized parameters defining the TD equilibria used in the paper: R 
and a are the torus major and minor radii, respectively, d is the depth 
of the torus center, L is the distance of the magnetic charges to the 
center of the torus, and zo and m are the depth and strength of the 
magnetic dipoles, respectively. 



sen considerable depth of the dipoles. Second, the radial profile 
of the current density is changed such that it no longer peaks 
at the surface but in the middle of the channel, and decreases 
monotonically towards zero at the surface. This removes unreal- 
istically steep gradients of the current density at the surface of 
the torus, including the magnetogram plane. The expressions for 
the equilibrium with this parabolic current density profile will be 
provided in a future publication (V. S. Titov, in preparation). 

We consider the modified TD equilibrium for several param- 
eter sets, four stable configurations and an unstable one. The pa- 
rameters defining each equilibrium are given in Table [T] All vari- 
ables are normalized by a characteristic length, taken here to be 
the apex height of the geometrical torus axis, R-d,hy the mag- 
netic field strength at the apex point, and by quantities derived 
thereof. 

Three of the four stable equilibria have identical geometric 
flux rope parameters and an almost identical left-handed, av- 
erage end-to-end twist of ^ 2. 1 tt in the center of the current 
channel (see Sect. |4] for the averaging procedure), but difl'erent 
properties of the ambient field. In the first two cases an HFT is 
present, staying close to the photosphere in one (Low_HFT) and 
reaching significantly into the volume in the other (High_HFT). 
A third case (No_HFT) has no HFT above the photospheric 
plane, hence, it has a simpler magnetic structure than the cases 
with an HFT. Due to the relatively strong toroidal field compo- 
nent chosen, these three configurations do not have bald patches. 
The fourth stable case (BP) has a left-handed average twist of 
^ I.Stt, close to the twist of the first three equilibria, but the 
minor radius of the torus is enlarged (and ^et correspondingly 
reduced) to introduce a bald patch in the resulting field. 

Figure [T] illustrates for the High_HFT and BP cases that 
the field strength in the magnetogram decreases rapidly with 
distance from the four flux concentrations. The distribution of 
a(x, y,z = 0) in the same figure shows that the field is highly 
nonlinear (see below for the occurrence of currents outside the 
torus in these plots). The transverse field near the section of the 
neutral line between the flux concentrations is strongly aligned 
with the neutral line for the BP configuration, as in a filament 
channel on the Sun. A zoom into this area is plotted in Figure [2l 
demonstrating by the inverse crossing of the neutral line in the 
central part of the magnetogram that a BP is indeed present. 

To obtain an unstable case, we increase the average twist to 
^ 2.7;: by reducing the minor torus radius. This configuration is 
unstable with respect to the helical kink instability. 

The analytical expressions in TD and their modifications de- 
scribed above represent only approximate equilibria. Moreover, 
the numerical determination of the current density introduces ad- 
ditional ("discretization") errors. Therefore, the individual mag- 
netic configurations were first relaxed to numerical equilibria. 
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(c) 



(d) 





Fig. 3. Selected field lines for the Low_HFT case; (a) reference, (b) extrapolated, (c) potential, and (d) linear field (with a = -0.2). 
In (a) and (b) the flux rope is depicted by two groups of field lines (in red and green), each one starting from six uniformly spaced 
points on a circle of radius 0.1, centered at the intersections of the magnetic axis of the flux rope (CFL) with the z = plane. Grey 
field lines are started from the z-axis and visualize the potential field above and below the current channel. In (c) and (d) the field 
lines are started from the same set of footpoints as in (a) and (b). Contour plots of Bz(x,y, 0) and its neutral line B^ix^y, 0) = 
(dot-dashes) are shown in the bottom plane of each panel. 



using the MHD code described in lTorok & Klieml (l2003h on uni- 
form Cartesian grids with a resolution A = 0.06 in all directions. 
This introduces only little change in the shape of the flux rope, 
but new layers of current density, not present in the analytical 
model, form in the separatrix and its vicinity, primarily below 
the current channel. These are visible in the a(x, y, 0) plots of 
Fig. [T] as narrow strips of forward and reverse current between 
the footpoints of the channel. 

From each numerically relaxed stable TD equilibrium, a data 
cube of size [-3.03, 3.03] x [-4.95, 4.95] x [-0.06, 4.44] was ex- 
tracted, containing the nonlinear magnetic field to be used as 
the test field for the extrapolation code (the somewhat larger 
cube used for the unstable case is given in Sect. 14.31) . We re- 
fer to these data cubes as the reference field, B^^^, for the given 
set of parameters. The reference fields were then reconstructed 
by magnetofrictional extrapolation, employing the same grid as 
in the MHD relaxation. The High_HFT reference field was also 
constructed and reconstructed using a uniform grid resolution 
of 0.03 in each direction. Since the achieved accuracy of recon- 
struction turned out to be comparable to the one obtained with 
resolution of 0.06, we selected the latter value for the present 
study. Although the nonlinear solution is known in the whole 
discretized volume, and specifically on its lateral and top bound- 
aries, the vector field in the bottom plane, B^^^(x,y,z = 0), is 
the only information used in the reconstruction of all test fields. 



Table 2. Compatibility of magnetogram and extrapolation 
method 



Case 



^flux ■ 



"W 



^force ' 



lo^ 



To^ 



High_HFT 


-0.99 


2.60 


0.93 


Low_HFT 


-1.83 


5.84 


1.24 


No_HFT 


2.28 


8.56 


1.66 


BP 


3.52 


29.71 


2.47 


Unstable 


-4.21 


90.42 


0.12 



Normalized flux, force, a nd torque imbalance (as defined in 
IWiegelmann et al.|[2006bh of the five magnetograms analyzed. 



because this is the challenge that the extrapolation codes have to 
face when they are confronted with observed magnetograms. 

Since the reference TD fields satisfy Eqs. ([T]) only to the ac- 
curacy of the finite-diff'erence scheme used in the relaxation, 
it is necessary to quantify the compatibility of their magne- 
tograms with the force-free hypothesis implicit in the extrapo- 
lation method. For this purpose, we use normalized integral av- 
erage s of flux, force, a nd tor que in the lower boundary as defined 
in Wiegelman n et al.l (|2006b), which vanish for perfect compat- 
ibility. These metrics, reported in Table [2l show that the flux 
imbalance is negligible and that the residual force and torque 
values remain small, even for the unstable case. 



G. Valori et al.: Testing magnetofrictional extrapolation with the Titov-DemouUn model of solar active regions 



3. Extrapolation method 

The implementation of magnetofrict ional extrapolation used i n 
this paper was described in detail in 'V alori et aP (l2QQ5l 120071) . 
so that we can restrict ourselves here to a short outline of the 
method and its recent improvements. Our formulation of the 
method consists in a pseudo-temporal evolution of an initial field 
subject to the equation 



dB 

^- = V X (v X 5) + CL^^iV • B), 
ot 



where 



V := Cy- 



JxB 

B^ ' 



J = VxB, 



(2) 
(3) 



and cy and cl are numerical factors fixed by suitable stability 
criteria. The initial field is usually chosen to be a potential field 
derived from the "line-of sight" magnetogram Bf^(x,y, 0), with 
the full vector magnetogram B^^^(x, y, 0) overwritten in the pho- 
tospheric layer. The rightmost te rm in Eq. Q is the V • 5 cleaner 
(see, e.g. jKeppens et al.l 120031] and references therein). By ap- 
plying the divergence operator, it is readily seen that Eq. ^ re- 
duces to a diff'usion equation for the scalar field V B. Since also 
the fi rst term on the right hand side of Eq. Q is of parabolic 
type (ICrai g & Sneydlll98 6), Lorentz forces and errors ofV-B 
diff'use out of the computation domain through the boundaries in 
the course of the extrapolation. This is best seen by computing 
the time derivative of the magnetic energy £^niag = f dV B^/2. 
Multiplying Eq. ^ with B and integrating by parts, it is found 
that 



(4) 



^.^mag = -CyJdV \JJ^ - Cl J dV (V • B)\ 

which implies that the magnetic energy decreases as long as 
there are Lorentz forces and magnetic sources in the domain. 
Here it has been assumed that the boundaries are perfectly 
conducting (as in ICraig & Sneydl 11986 ) and that the field is 
solenoidal at them, so that surface integrals do not contribute. 
During the pseudo-temporal evolution, the decrease of the per- 
pendicular current density J^ reduces also the relaxation ve- 
locity V [Eq. ([3])]. If a static state is reached, then the field is 
divergence-free and v = implies that it is also force-free. By 
construction, such a static solution extends smoothly from the 
vector magnetogram in the bottom boundary and represents the 
extrapolated nonlinear force-free field. We note that by overwrit- 
ing the vector magnetogram in [z = 0}, we actually introduce a 
finite jump in the field close to the photosphere, yielding Lorentz 
forces and errors in the solenoidal property; consequently, sur- 
face integral terms in Eq. Q are important initially. This results 
in an initial rise of the energy not described by Eq. ^. A discus- 
sion of this eff'ect will be given in a future paper. 

Equatio n Q can be form ally derived from the MHD equa- 
tions (as in'Yanget al.l ll986b by assuming that viscosity dom- 
inates the evolution (hence the characterization of the method 
as "magnetofrictional"). However, since such a strong viscosity 
has no immediate physical justification for the nearly collision- 
less coronal plasma, one can equivalently consider Eq. (O to be 
a numerical tool that finds a force-free magnetic field for given 
boundary conditions, with only the initial potential field (before 
the vector magnetogram is overwritten) and the final nonlinear 
force-free field having a well defined physical meaning. 

Our implementation of the magnetofrictional code employs 
a fourth-order central-diff'erence scheme for the space discretiza- 
tion. The pseudo-time discretization is based on a gene ral MHD- 
stability prescription (see details in I Valori et al.ll2007h . In order 
to reduce the execution time, the time discretization was recently 



compl emented by a Runge-Kutta-Chebyshev acceleration tech- 
nique ('Alexiades eraDll996l:lEvie et al.l l200'^. The detailed de- 
scription and testing of this upgraded pseudo-time discretization 
is beyond the scope of the present paper and will be provided in 
a future publication. 

The field in the photospheric boundary is fixed to keep the in- 
put magnetogram values throughout the whole extrapolation. On 
the side and top boundaries of the extrapolation volume, which is 
identical to the volume of the reference field given in Section O 
the normal component is derived from the inner field values so 
to have vanishing V • Bin the boundary. The transverse compo- 
nent is given by a fourth-order polynomial extrapolation of the 
inner field values o nto the boundaries. These boundary condi- 
tions were tested in I Valori et al.l (l2007h . where more detail can 
be found. Finally, all extrapolations presented here start from 
the potential fi eld obtained by Schmidt's method (ISchmidtl 19641: 
ISakuraill 19891) . except for a potential-field and a linear-field case 
discussed in Sect. 14.21 for which the method by lSeehafe3 (Il978h 
was used. 



4. Analysis of the extrapolation results 

The reconstructed fields, B{x, j, z), are compared to their respec- 
tive reference field, 5^^^(x,j,z), in an "analysis domain" given 
by clipping 20 grid layers from the top and lateral boundaries of 
the extrapolation volume. This is done in recognition of the fact 
that the magnetograms have considerable margins of current- 
free field around the flux rope footpoints. The currents remain 
relatively small in the corresponding parts of the extrapolation 
volume, where the relaxation is slowed and, therefore, the qual- 
ity of reconstruction is reduced. Sim ilar practic e, with larger 
margins, was adopted in lSchrijver et al. (2006) and Metcalf et al.l 
(2008). The quality of the reconstructed equilibria is quantified 
in the following in several ways. 

The degree of force-freeness is measured by the current- 
weighted, average sine-angle between current and magnetic 
field. 



crj = 



1.1 Ji ' 



where cr, = 



BiJi 



(5) 



and the sums run over all grid nodes / in the analysis domain. 
This sam e quantity is sometimes referred to as CWsin (e.g., 
ISchrijver et al. 2008a). The solenoidal property is quantified, 
following i Wheatla nd et al. (2000), by the average over the grid 
nodes {\fi\) of the fractional flux 



fi = 



UdS-B, ' 



(6) 



through the surface d*V of a small volume *V including the 
node /. For uniform grid spacing A, this reduces to fi = A(V • 
Bi)/(6\Bi\) and < \fi\ ) = 2/ \fi\/N, where A^ is the number of grid 
points. The influence of the residual spurious magnetic charges 
implied hyV-B^O can also be quantified by the non- solenoidal 
part of the magnetic field, ^ns- Using Gauss' theorem, we have 



^n 



A^ 
47: 



X(v-B)-- 



X - Xi\- 



(7) 



which is evaluated on an auxiliary grid displaced by A/2 in each 
direction. All three qu antities shoii l d be a s small as possible. We 
have demonstrated in IValori et aP (l2007l) that our code reduces 
the value of V • 5 to machine precision in some cases; however. 
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due to the parabolic nature of Eq. (|2]), this may require long ex- 
ecution times. For the purposes of the present work, we termi- 
nated the relaxation when the force-freeness metric crj tended 
to reach a plateau of sufficiently small value (of order 10"^) and 
{\fi\) was of the same order of magnitude as the corresponding 
value for the initial potential field or smaller. 

Since the reference field B"^^^ is known, we can perform a 
direct comparison between B^^^ and the reconstructed field B 
obtained from t he magnetograni. Sev eral comparison metrics 
were defined in ISchrijver et al.l (l2006h . out of which we adopt 
the mean vector error 



]_y\Bf-Bi\ 



\B 



ref\ 



(8) 



The smaller £'m, the better the reconstruction. As for previous 
reconstructions of te st fields ( Schrijver et al. 2006; Valori et al. 
l2007l: iMetcalf et al.l l2008), this metric is the most sensitive one 
to reconstruction errors of th e TD field. The furthe r metrics of 
direct comparison defined in ISchrijver et al.l (l2006h . Cyec, C'cs, 
£^N, yield even smaller deviations from the ideal value (1 or 0) 
than Em for all cases analyzed here. 

To judge geometry and topology, we first determine the apex 
height of the flux rope's magnetic axis, also referred to as the 
"central field line", CFL, in the following. (The CFL diff'ers 
slightly from the geometrical axis of the torus, since the exter- 
nal potential field drops with height and, hence, contributes dif- 
ferently to the bottom and top parts of the flux rope.) Next we 
check for the presence of an HFT and its apex height. Both apex 
heights are given by inversion points of the poloidal component 
of the field at the line- symmetric z axis. This is well approx- 
imated by Bxifi, 0, z) because the flux rope writhes only rather 
weakly out of the plane [x = 0} in the equilibria considered in 
this paper. Finally, we check whether the BP in the fourth stable 
test field is recovered at the correct section of the neutral line. 

The TD eq uilibrium can be susceptible to the helical 
kink instability (T5roket al .1120041) and to the torus instability 
(iKliem & Torokll2006b . The parameters controlling these insta- 
bilities should be reliably and accurately reproduced by the ex- 
trapolation. The main parameter which regulates the stability 
with respect to the helical kink mode is the field line twist. This 
quantity, O = lB(^(p)/(pB^(p)), is defined with respect to the 
CFL of length / using local cylindrical coordinates (p, 0, f ). We 
use its average in a circular cross section of normalized radius 
0.3 centered at the apex of the CFL. Within this radius, the av- 
erage twist is ( O ) ^ 2.1 TT for three of the four stable equilibria 
(High_HFT, Low_HFT, and No_HFT) and ^ L8 ;r for the fourth 
stable case (BP). For the unstable case discussed in Sect. 14.31 
the average twist in p < 0.3 is ^ 2.7 n. The threshold of the 
torus instability is given in terms of the "decay index" of the ex- 
ternal poloidal field, n = -RdlnB^^ldR, where R is the major 
torus radius. F or n > Hex, with the threshold lying in the range 
^cr ~ 3/2 ... 2 (JKliem & Torok 2006), an expansion of the cur- 
rent ring out of equilibrium cannot be balanced by the external 
poloidal field. Computing the index n requires the knowledge of 
the external poloidal field, which is not immediately available 
for our numerical solutions, since it cannot be easily separated 
from the poloidal field produced by the ring current. Therefore, 
we compute the decay index of the total poloidal field at the z 
axis, given to a good approximation by 5^(0, 0, z). 



n = i-z 



dlnBMO^z) 
~Jz 



(9) 



1.55 


0.73 


1.18 


0.80 


6.04 


6.90 


10.28 


10.05 


0.044 


0.036 


0.033 


0.054 


1.94% 


0.54% 






4.23% 


0.63% 


2.70% 


0.38% 


2.14% 


1.10% 


2.69% 


0.44% 


3.92% 


3.89% 


6.61% 


0.51% 


0.22% 


0.37% 


0.74% 


0.03% 


4.55% 


4.43% 


4.83% 


2.80% 



The average is taken in a small interval in z, which is located just 
above the current channel in all equilibria considered here. The 



Table 3. Reconstruction accuracy of the stable TD equilibria 



Figure of merit High_HFT Low-HFT No-HFT BP 

o-j X 10' 
(l/DxlO^ 

Em 

HFT apex 

CFL apex 

(cD> 

h 

^mag 
-'^mag 

Figures of merit in the analysis domain as defined in Sect. ID with 
smaller values indicating higher accuracy. Deviations from the 
corresponding values in the reference fields are given in percent. 



quantity h is not suitable to determine the onset of the torus insta- 
bility, but it serves the purpose of quantifying how well the sta- 
bilizing potential field above the current channel is reproduced 
by the extrapolation. 

The magnetic energy, £'mag, and relative magnetic helicity, 
^mag, should also be accuratel y reproduced by the reconstructed 
field. We use the formulae by iDeVord (l2000i) to obtain an ap- 
proximation of ^mag for both the reference and the reconstructed 
fields. 

In order to reliably compare the eff'ects that topological dif- 
ferences have on the extrapolation quality, all computations pre- 
sented in this paper for stable TD equilibria were performed us- 
ing identical numerical parameters, rather than being individu- 
ally optimized for best reconstruction. 

4.1. Stable equilibria 

We consider the reconstruction of four stable TD equlibria, two 
of which contain an HFT (of diff'erent heigh t) and one conta ins 
a BP (Table [B, using the potential field bv lSchmidtl (Il964l) as 
initial condition. The main results are summarized in Table O 
where measures of the reconstruction quality of the four refer- 
ence fields are listed (see also the comprehensive compilation of 
the various metrics for all our extrapolations in Table [5]) . The ta- 
ble shows that the errors are very small and of comparable mag- 
nitude for all cases (within fractions of a percent up to a few 
percent). In particular, we point out that the mean vector error 
Em of all four reconstructed fields is 0.054 or smaller, a result 
far better than the accuracy achieved in any previous similar re- 
construction stu dy, i.e., in those that used only the vector magne- 
togram as input (Schriiver et al. 2006: Wiegelman n et al.ll2006al: 
Valori et al._ _2007: Metcalf et al. 2008), although a smaller mar- 
gin between the analysis and extrapolation volumes is used here. 
This mean vector error corresponds to a mean angular deviation 
between B and B^^^ of at most 3 degrees (if the field strengths 
agree on all grid nodes) and to a mean deviation of field strength 
of at most 5% (if all directions agree) — values lower than the 
noise level for most pixels in observed magnetograms. It is also 
worth noting that the reconstruction accuracy is now compa- 
rable to that achieved previously in "extrapolations" that have 
provided the reference field vector on all six boundaries of the 
reconstruction volume. The accuracy reached here is clearly su- 
perio r to all but on e of th ose results reported in Tables 1 and 
2 in Schrijver et al.l (l2006l) and is approachi ng the best one in 
this paper, as we l l as th e results in Amari etaP (|2006) and 
IWiegelmann et al.l (l2006ah . The field line plots in Fig. [3] show 
that such small errors have practically no visible influence on 
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Fig. 4. Plane averages of the relative solenoidal error l^nsl/l^l in 
the No_HFT case, as a function of position. 



the reconstruction, which is virtually indistinguishable from the 
reference field. 

The best reconstruction is obtained for the case with a bald 
patch, while the least accurate one is the No_HFT case, where 
neither an HFT nor a BP are present. There is apparently no re- 
lation between the quality of the reconstruction (Table O and 
the compatibility of the corresponding magnetogram with the 
force-free constraints (Table [2]). This shows that the level of in- 
consistencies in the vector magnetograms of our stable reference 
fields is sufficiently small to have essentially no influence on the 
extrapolation quality. 

Turning to a more detailed evaluation, we first consider the 
consistency measures of the reconstructed fields given in the up- 
per part of Table [3] On average, the angle between the magnetic 
field and the current density stays below 1 degree, implying that 
the reconstructed fields have reached a high degree of force- 
freeness. The errors in satisfying V • 5 = as measured by ( |/i| ) 
are also very small, well below those of the potential field used 
as initial condition (which is nominally divergence-free, except 
for discretization and truncation errors). Similarly, l^nsl/l^l is 
small. Let us consider this quantity for the No_HFT case which 
has the highest value of < |/i| ). Since the directions of ^ns(-'^) 
are essentially random, a very conservative estimate of its influ- 
ence is ( l^nsl/l^l ), which stays below 6%, with a grid average 
of 0.7%. Averaging the relative error l^nsl/l^l in the planes of 
constant x, j, and z, yields the profiles plotted in Fig. IH These 
show that the deviations from the solenoidal condition originate 
primarily at the boundaries of the box, actually being smallest 
in the current channel in its interior. Overall, the values of ct/, 
( ly^l ), and B^^I\B\ demonstrate that the reconstructed stable TD 
equilibria possess a high degree of consistency with the assump- 
tions of the extrapolation method. 

The geometry of the reference field is perfectly recovered, 
as expected from the very low value of Em and illustrated in 
Fig. [3] for the case Low_HFT. The flux rope, as the entire ex- 
trapolated field, reproduces the z-axis line symmetry of the TD 





Fig. 5. Isosurface of the reconstructed current channel at 3 1 % of 
peak current density for the cases High_HFT (top) and No_HFT 
(bottom). Field lines crossing the z axis are shown in grey if be- 
low the HFT and in green if above the HFT. 



equilibrium exactly. Therefore, the flux rope location is correctly 
represented by the apex height of the CFL, found to match the 
height in the reference field nearly exactly for all stable equilib- 
ria. Since Em includes a normalization to the local field value 
on the grid, it gives the weak-field areas in the outer parts of 
the reconstruction volume the same weight as the strong-field 
areas near the bottom and in the flux rope. From the very small 
value of Em, we can therefore be sure that the high overlying 
field loops are reconstructed with essentially the same accuracy 
as the flux rope, and this is apparent from Fig. [3] as well. The 
comparison of such loops with soft X-ray or EUV loops has 
proven insatisfactory in rece nt nonlinear extrapo lations of ob- 
served vector magnetograms (IPeRo sa et al.ll2009l) . although the 
observed loops very likely outline force-free or even current- free 
field line bundles. The results given here demonstrate that this is 
not due to limitations of the extrapolation scheme(s). 

The HFT is recovered for both reference fields containg one, 
with high accuracy of the HFT apex height. This is illustrated for 
the High_HFT case in Fig. [51 which also includes the No_HFT 
case for comparison. The bald patch in the BP case is recovered 
as well, occupying a similar section of the neutral line as in the 
reference field, see Fig. [2l These results demonstrate that our 
magnetofrictional code is able to reconstruct the geometry and 
topology of TD-like equilibria. 

Considering the parameters that characterize the relevant in- 
stabilities, we find that the relatively high twist of ^ (1.8-2.1);: 
is recovered with very high accuracy, the error s being below 3% 
in all four cases. Thus, as in lValori et al.l (2005), it is confirmed 
that the magnetofrictional method can reconstruct field lines that 
wind around a flux rope axis more than one time. In this regard 
the nonlinear extrapolation difl'ers in principle from linear ex- 
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Fig. 6. Poloidal field at the z axis, Bj^ifd, 0, z), in the No_HFT case, 
for the reference (red), extrapolated (green), and TD_ex poten- 
tial fields (blue). Left: entire z axis extension; right: upper part 
(above the current channel) in logarithmic scale. Note that the 
off'set of the blue curve results from normalizing the external 
field in the same manner as the full TD field, i.e., by the value 
|5(0, 0, 1)1, even though it does not contain a flux rope. 



trapolations, which limit the number of field line turns in the 
computed field to about 1/2 for realistic magnetogram sizes, due 
to the upper limit on the force-free paramete r inherent in linea r 
fields (Kliem et al., in preparation; see also 'Leka e t al.ll20Q5b . 
This means that estimates of twist in fields extrapolated with the 
magnetofrictional code are by far more reliable than estimates 
using a linear code. 

In order to assess the reliability of the extrapolation with re- 
gard to the torus instability, we consider the approximate decay 
index of the total poloidal field, h, given by Eq. ©. The No_HFT 
case is illustrated in Fig.[6]to show that the reference and recon- 
structed profiles of Bx{^, 0, z) agree nearly exactly in the whole 
range of z values. The same is true, of course, for the slope h, 
especially in the relevant region immediately above the current 
channel, as the logarithmic plot on the right-hand side of the fig- 
ure shows. The very good agreement of the field profile at the 
whole z axis implies that both parts of the poloidal field must be 
well recovered individually and confirms the use of h as proxy 
for the reconstruction quality of n. Figure [6] suggests the height 
range above the current channel not much exceeding z = 2 for 
the computation of a representative value for h. We will use the 
range z e [1.80, 2.16]. The values thus obtained diff'er by only 
0.5-7% from the true ones (Table [3]), a precision better than that 
of the current knowledge of the torus instability threshold. The 
comparison with the corresponding profile of the external field 
components in the TD equilibrium (blue curves) illustrates that 
the slopes of the field, n and h, indeed diff'er considerably in the 
height range used for the h computation, so that the absolute val- 
ues of this parameter, included in Table[5]for completeness, have 
no direct bearing on the torus instability. Since the close agree- 
ment of h for the reconstructed and reference fields implies that 
n is recovered with similar accuracy, we conclude that the key 
parameters which determine the onset of the two relevant insta- 
bilities are reliably reproduced. 

Similarly, the magnetic energy and relative magnetic helicity 
are obtained with high accuracy; their errors remain below 1 and 
5 percent, respectively, for all four stable reference fields. 

Figure [3] also shows that the potential and linear extrapola- 
tions entirely miss the flux rope, typically leading to an arcade- 
type structure. Despite the failure of the potential field to re- 
produce any of the salient features of the reference field, still 
some of the comparison metrics score decently (e.g., Em = 0.28; 
see Table [5] for more detail). This is due to the fact that a large 



fraction of the analysis volume is occupied by potential field. 
The energy of the potential field reconstruction, being the min- 
imum for a given normal field distribution at the boundaries, is 
of course lower than the energy of the reference field. However, 
it still accounts for a large part of it, about 80% in the selected 
analysis volume. The linear reconstruction suff'ers from similar 
limitations but, since it fills the entire volume with currents, it is 
aff'ected by even larger errors than the potential field. The value 
a = -0.2 adopted for the computation of the linear field is the 
one that gives the best r.m.s. matching of the tran sverse magne- 
togram components, often referred to as o^best (see lPevtsov et al.l 
1995). 

Finally, w e briefly comment on the re construction of the TD 
equilibrium in lWiegelmann et aP (l2006ah . The original TD equi- 
librium with a line current was considered, and a difl'erent ex- 
trapolation scheme, the optimization method, was used. The pa- 
per includes a reconstruction based only on the information from 
the bottom boundary as "Case II." The corresponding field line 
plot in Fig. 2 of that article (done by one of the authors of the 
present paper, GV) indicates that the extrapolated field contains 
two partially twisted flux tubes rather than one. That plot is ac- 
tually misleading because the appearance of two separate flux 
tubes has been obtained by calculating field lines starting around 
the nominal footpoint positions of the geometrical torus axis; 
however, as described above, the footpoints of the flux rope's 
magnetic axis do not coincide with those of the geometrical torus 
axis. A proper field line plot (not included here) is obtained by 
placing the field line start points around the magnetic axis, and it 
shows a single, higher and slightly writhed flux rope, similar to 
the cases analyzed here. Due to the many difl'erences between the 
employed test fields and reconstruction volu mes, we do not fur- 
ther pursue the comparison with the results in lWiegelmann et aP 
(2006a). Those are far less accurate than the ones given here in 
the case that only the vector magnetogram is used as input. 

4.2. Dependence on the initial field 

In this section we describe the link between the potential or lin- 
ear force-free field that is used as initial condition for the ex- 
trapolation and the quality of the obtained reconstruction. A key 
consideration here is the balance between the hoop force and 
the Lorentz force provided by the external poloidal field. We 
consider the Low_HFT case in combination with four diff'erent 
initial fields: tw o p otential fields c omputed by the methods of 
ISchmidtl (11964 and lSeehafe3 (Il978l) . respectively, a linear force- 
free field computed by the Seehafer method, and the external 
potential field of the TD equilibrium (generated by the magnetic 
charges and dipoles in the modified version of the equilibrium 
used in this paper). The latter, referred to as TD_ex in the fol- 
lowing, is, of course, not available for observed magnetograms. 
Figure [7] shows the poloidal component at the z axis for the 
reference field (red line), for the four fields used as initial condi- 
tion (blue lines), and for the corresponding four reconstructions 
(green lines). The Schmidt and Seehafer potential fields are, of 
course, computed using the same Bz(x,y, 0) as boundary condi- 
tion, but they diff'er in their assumptions about the field outside 
the magnetogram area. While Schmidt's method assumes the 
field to vanish outside (and therefore requires flux balance within 
the magnetogram), Seehafer' s method uses a periodic repetition 
of an extended magnetogram obtained by mirroring the original 
magnetogram about one corner (in this way ensuring the flux 
balance). The resulting poloidal component of the Schmidt po- 
tential field approximates the external poloidal field of the TD 
equilibrium rather well, being larger by only a factor 1.06 at the 
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Fig. 7. Poloidal field at the z axis, 5^^(0,0, z), in the Low_HFT 
case, for the reference (red solid line), initial (blue), and cor- 
responding extrapolated fields (green), for diff'erent type of the 
initial field: potential Schmidt (short dashes), potential Seehafer 
(dot-dash), linear alphabest Seehafer (triple dots-dash), and po- 
tential TD_ex (long dashes). 

Table 4. Reconstruction accuracy of the Low_HFT equilibrium 
for diff'erent initial conditions 
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nominal position of the CFL apex, z= 1.14. Hence, the flux rope 
to be built by the extrapolation is provided with almost the cor- 
rect external poloidal field, so that the force balance is indeed 
found very close to the correct position (Table |4]). 

The Seehafer potential and linear fields are weaker at the 
nominal CFL apex position, with a ratio of the poloidal compo- 
nents to ^ep of the TD_ex field of 0.9 and 0.48, respectively, so 
that the flux rope finds equilibrium at greater heights. However, 
the apex heights difl'er by only 14% and 18%, respectively, 
which shows that the initial field does have an influence on the 
extrapolation, but also that the magnetofrictional relaxation can 
tolerate a substantial mismatch between the initial field and the 
true external field while still producing relatively small errors. 
The information flowing from the magnetogram into the extrap- 
olation volume can largely, but not completely, modify the prop- 
erties of the initial field, although open boundaries are imple- 
mented at the top and sides of the box, which allow the field to 
change in response to changes in the interior. This has to do with 



the fact that large parts of the volume are current free, both in 
the initial and in the final configuration. Hence, currents must 
be built up and subsequently be removed, in order to modify the 
initial field in the large outer parts of the box. This is a slow 
process, eventually hindered by numerical diff'usion. 

Apart from the CFL apex height and the related decay in- 
dex h, the figures of merit of the four extrapolations in Table IH 
show excellent values. The tendency for the Seehafer potential 
and linear fields to score worse is visible but rather weak, indi- 
cating that the initial field has a weaker influence on quantities 
like flux rope twist and total energy and helicity. These depend 
primarily on the information contained in the magnetogram and 
only to a smaller extent on the flux rope height. It is also clear 
that the figures of merit for the Schmidt potential field mostly 
beat those of TD_ex because the former is more consistent with 
the magnetogram, which includes the contribution by the current 
channel. 

All four extrapolations recover the very low lying HFT apex 
of the test field, essentially at the correct height (Table [5]) . 

These comparisons clearly suggest to use the Schmidt po- 
tential field as initial condition for flux-balanced magnetograms 
of isolated active regions (two conditions that go together in 
tendency). A flux imbalance can often be reduced by embed- 
ding the vector magnetogram in a larger magnetogram, which 
is typically of lower resolution and (with current instrumenta- 
tion) only a line-of-sight magnetogram. The transverse compo- 
nents must then be modeled, for example by assuming a poten- 
tial field, and this can strongly influence the outcome of the ex- 
trapolation. Alternatively, the use of the imbalanced vector mag- 
netogram joint with open boundary conditions that allow for flux 
and currents to leave the numerical box through the top and side 
boundaries, as in ou r code, can lead to better extrapolations (see 
iDeRosa et al.l 12009 and Fuhrmann et al. 2010 [in preparation] 
for a discussion). Finding the best trade-off* between these op- 
tions will require further study. 

4.3. Unstable equilibrium 

As mentioned in the Introduction, the parameters defining the 
TD equilibrium can be chosen such that the configuration is un- 
stable with respect to the helical kink or to the torus instabil- 
ity. Such equilibria off'er the opportunity to test the extrapolation 
code in the reconstruction of unstable force-free configurations, 
addressing the question which signatures of the unstable situa- 
tion are produced. We compare the extrapolation of the magne- 
togram provided by an unstable TD equilibrium with the equi- 
librium itself and with the post-eruption configuration obtained 
by a numerical MHD evolution of the perturbed equlibrium with 
the magnetogram kept fixed. 

The configura tion we con sider is the same as the eruptive 
case in iTorok & Klieml (12005), except that a smaller initial twist 
of 2.7;: is used (note that the twist average ove r the whole cross 
section of the current channel, as calculated by Torok & Klieml 
is ^ 47: for this configuration). The MHD simulation is per- 
formed in a Cartesian box of size [-20, 20] x [-20, 20] x [0, 40], 
which consists of the uniform grid used for the extrapolation 
in its central part and a streched grid outside the extrapolation 
domain. The smaller twist allows us to relax the configuration, 
including the magnetogram plane, for a few Alfven times, reduc- 
ing the initial spurious Lorentz force densities, before the helical 
displacement due to the developing kink instability becomes vis- 
ible. 

After this initial relaxation, we fix the magnetic field in 
{z = 0} for the whole subsequent MHD evolution. We refer to 
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Fig. 8. Selected field lines for the reference fields of the unstable 
case at the initial (top) and final (bottom) stages of the MHD 
evolution (see Fig.[3]for the color coding of the field lines). 



this partially relaxed, unstable configuration as the initial refer- 
ence field for the extrapolation of the unstable case; a field line 
plot is provided in Fig.[8l A small, transient, upward directed ve- 
locity perturbation is then applied at the flux rope apex to ensure 
that the rope kinks upwards. Once the kink instability has lifted 
the rope to a height where the field decreases sufficiently fast, 
the torus instability sets in. The rope is then additionally accel- 
erated by the torus instability and erupts fully. In the wake of the 
rising flux rope, a vertical current sheet is formed. Reconnection 
in this current sheet cuts the legs of the flux rope at low heights 
(whereas the bulk of the rope still erupts), and the legs then con- 
nect to the sunspots. The upper part of the reconnected flux rope 
reaches the closed top boundary at ~ 90 Alfven times and is sub- 
sequently compressed there. In the remaining time of the simula- 
tion, which is terminated after 750 Alfven times, changes in the 
lower part of the box occur mainly in the cusp- shaped arcade be- 
low the reconnection region, whereas the reconnected flux rope 
legs remain practically unchanged. 

In order to be compatible with fixing the magnetogram, we 
set the plasma velocities to zero in the grid layers [z = 0} and 
{z = A}. Therefore, although the simulation is run for a long 
time, the Lorentz force densities do not relax completely, so that 
the post-eruption configuration, although largely relaxed, is not 
perfectly force-free in the lower layers of the box. Nevertheless, 
the field line connectivities do not change significantly anymore, 
so that a qualitative comparison with the force-free extrapolation 
is possible. We refer to this post-eruption configuration as the 
final reference field for the extrapolation of the unstable case 
(see again Fig. [8]). 

The extrapolation of the unstable case is obtained in a sim- 
ilar way as described in Sects. |4J1 and 14. 2 1 except that here the 
extrapolation grid discretizes a larger volume [-3.00,3.00] x 
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Fig. 9. Evolution of crj (solid line) and £^mag (dashes) during the 
extrapolation runs of the cases BP (left) and Unstable (right). 
The energy is normalised to its value at the first iteration. 



[-5.04,5.04] X [-0.06,6.12] with the same uniform resolution 
A = 0.06. The extrapolation employs the z = magnetogram of 
the initial/final reference field and starts from the Schmidt po- 
tential field in [z > 0}. Table [2] verifies that the magnetogram 
includes residual forces, but it has a lower value of the applied 
torque than the stable cases. According to the metrics in the 
table, the unstable near-equilibrium TD configuration obtained 
from the initial partial MHD relaxation (before the instabilities 
become dominant) provides a magnetogram that is indeed ac- 
ceptably force free. 

In all stable cases analyzed in the previous sections, the mag- 
netic energy and crj evolve monotonically during the whole ex- 
trapolation run (£^mag increases, crj decreases), until nearly con- 
stant values are reached in ~ 7000 iterations. The magnetic field 
hardly changes after this point, even if the magnetofrictional re- 
laxation is continued for a very large number of iterations. An 
example of this behavior is given in Fig. [9] for the BP case. The 
extrapolation of the unstable case is significantly diff'erent. As 
Fig. [9] shows, the energy (resp., crj) reaches a temporary maxi- 
mum (resp., minimum) in about the same number of iterations, 
but then the field changes and rapidly reaches an almost flat 
plateau with a lower energy and higher crj. This plateau extends 
over -2x10^ iterations, followed by a very gentle "bump" ex- 
tending over -5x10^ iterations, and then by a second plateau. 
In the whole evolution from the beginning of the first plateau, the 
field line plots show that there is hardly any change in the mag- 
netic configuration, so that we terminated the run after 9 x 10^ 
iterations. From the extrapolation run we select a maximum- 
energy and a minimum- energy reconstruction, respectively cor- 
responding to the maximum value of the energy at 6700 itera- 
tions and to the almost static state reached at the end of the run 
(referred to as "MF Schmidt max" and "MF Schmidt min" in 
TableO. 

The maximum-energy reconstruction is a relatively force- 
free {(Tj = 0. 13) and divergence-free field, albeit the correspond- 
ing metrics reach only values that are about one order of magni- 
tude worse than for the stable cases (see Table O. Consequently, 
the r.h.s. of Eq. (|4]) does not vanish, although ^?£^mag = 0, which 
illustrates the limitations of the applicability of Eq. Q discussed 
in Sect. [3] This state is clearly a transitory one, as is the initial 
reference field. The two are determined by a transition between 
processes of partially diff'erent character, so that one can expect 
them to be neither very similar nor completely diff'erent. In the 
MHD run, the transition occurs between the relaxation of the 
initial Lorentz forces in a field that already includes a current 
channel and the dominant evolution of the instabilities, which 
release part of the free magnetic energy contained in the cur- 
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Fig. 10. Selected field lines for the extrapolation of the unstable case at diff'erent iterations. From left to right: the initial potential 
field at it=l, the maximum-energy configuration at 6700 iterations, and the final low-energy configuration at 9 x 10^ iterations. See 
Fig. [3] for the color coding of the field lines. 



rent channel. In the magnetofrictional run, the transition occurs 
between a phase of buildup of the current channel and the subse- 
quent relaxation to a minimum-energy state, which destroys the 
channel. The maximum-energy field contains a flux rope with an 
HFT underneath, but this rope is not as compact as in the initial 
reference field (compare the top panel of Fig. [8] with the middle 
panel in Fig. [T0|. Both the HFT and CFL reach about half the 
corresponding heights of the initial reference field. 

We recall that the extrapolation starts from a potential field 
that has obviously no flux rope (as shown in the left panel of 
Fig. [TOl) , and therefore the flux rope in the maximum-energy 
reconstruction is formed by the extrapolation process. The ex- 
trapolation does not succeed to form the flux rope completely 
because the underlying TD configuration is in the unstable do- 
main, so that the Lorentz forces tend to propagate the informa- 
tion about the current channel from the magnetogram into the 
volume and tend to move the system away from that configura- 
tion at the same time. The partial buildup of the current channel 
is also reflected by the fact that the maximum-energy field has 
reached 77% of the energy in the initial reference field. In light 
of this number, the very close agreement of the magnetic helic- 
ities must be considered a coincidence or a consequence of the 
approximate nature of its calculation (using the expressions in 
lDeVore,.200 0). 

Conversely, the minimum-energy solution has no flux rope: 
field lines started at the original location of the flux rope connect 
to the sunspots (rather than to the flux rope footpoints in the TD 
and initial reference fields), in a manner largely similar to the 
final reference field (compare the bottom panel in Figs. [8] with 
the right panel in[TO]). The magnetic energy has a smaller relative 
error (of 16%) than the maximum-energy reconstruction, while 
Em = 0.37 is a rather poor match. Although both the MHD and 
magnetofrictional runs eventually proceed to a low-energy state, 
the closed box used in the MHD run prevents the system from 
relaxing as deep as in the magnetofrictional run. Furthermore, 
as argued above, we cannot expect the reconstruction of the final 
reference field to be accurate, since this field presents relatively 



strong forces close to the magnetogram due to the MHD eruption 
process. 

5. Conclusions 

This investigation demonstrates that a force-free model of solar 
active regions, which includes relevant structural (i.e., topolog- 
ical or geometrical) features and is free of noise, can be recon- 
structed to a very high accuracy based only on its vector magne- 
togram. The reconstructed configuration consists of a flux rope 
with a non-neutralized current channel in its core. A bald patch 
separatrix surface or two quasi- separatrix layers combined into 
a hyperbolic flux tube form the interface between the flux rope 
and the ambient potential field. These elements are regarded to 
be generic for many active regions, so that the model realistically 
captures their structural properties, in particular for essentially 
bipolar flux distributions. 

Previous reconstructions of nonlinear force-free test fields 
considered the simpler configurations of an approximately 
current-n eutralized flux rope (Valori et al. 2005) and of a sheared 
arcade ( Schriiver et"aD 120061: [Valori et alj 120071) . A force-free 
model of a filament channel containing a hollow current chan- 
nel could also be reconstructed, albeit at significantly lower 
accuracy, which was likely due to the combination with a 
nois y line-of- sight ma gnetogram in the exterior of the chan- 
nel ( Metcalf et al.l l2008). Thus, the reconstruction capability has 
now been demonstrated, mostly at excellent accuracy, for all ba- 
sic types of force-free equilibria currently considered in coronal 
physics, with the one treated here providing the highest degree of 
structural co mplexity. This confirms and substantia tes previous 
conclusions (Metc alf et aD2008l:lDeRosa et al.|[200 9) that insuf- 
ficient accuracy of nonlinear force-free extrapolation originates 
from the violation of the force-free assumption by the input vec- 
tor magnetogram, not from insuflftcient intrinsic accuracy of the 
extrapolation codes. 

The accuracy achieved in the reconstruction of stable modi- 
fied TD equilibria clearly permits to discriminate between con- 
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Table 5. Extrapolations results for all TD equilibria 
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Reference 


0.0635 


1.137 


-2.128 


2.435 


1.0000 


0.3789 


10.087 


3.197 


4.303 


MF Schmidt 


0.0639 


1.130 


-2.121 


2.340 


0.9645 


0.7333 


10.050 


3.055 


6.899 


MF Pot Seehafer 


0.0624 


1.295 


-2.176 


2.008 


0.9355 


4.1202 


10.319 


3.241 


16.18 


MF Lin Seehafer 


0.0622 


1.341 


-2.192 


1.537 


0.9296 


4.5231 


10.315 


3.416 


14.47 


MF TD_ex 


0.0629 


1.183 


-2.183 


2.494 


0.9284 


1.1770 


10.237 


2.954 


11.98 
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1.928 
2.512 
4.130 


0.7246 
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0.6665 




7.958 
8.143 
8.267 






2.880 


277.9 
1.247 
1.225 


No_HFT 


Reference 
MF Schmidt 


1.131 
1.101 


-2.111 
-2.097 


2.350 
2.195 


1.0000 
0.9670 


0.2951 
1.1756 


10.047 
9.972 


3.105 
2.955 


4.254 
10.28 


BP 


Reference 
MF Schmidt 


1.389 
1.385 


-1.806 
-1.814 


-0.6210 
-0.6242 


1.0000 
0.9460 


0.0538 
0.8003 


7.489 
7.487 


3.610 
3.509 


7.162 
10.05 



Unstable 


Reference ini 


0.1705 


1.008 


-2.853 


2.916 


1.0000 


3.4610 


3.286 


24.21 


3.841 


Reference fin 








2.410 


1.0000 


22.116 


2.205 


23.69 


16.86 


MF Schmidt max 


0.0984 


0.4095 


-4.115 


1.992 


0.7693 


13.140 


2.540 


25.35 


94.96 


MF Schmidt min 










0.6266 


27.396 
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23.08 
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Reconstruction measures (absolute values) for all extrapolations presented in this paper, belonging to the cases No_HFT, Low_HFT, High_HFT, 
BP, and Unstable, supplemented by similar analysis for the diff'erent initial fields used in the Low_HFT case. The leftmost column characterizes 
the analyzed field using the following naming convention: 'Reference' identifies the model field to be reconstructed by the extrapolation; 'MF' 
stands for nonlinear magnetofrictional extrapolation, 'Pot' and 'Lin' for potential and linear extrapolation, respectively; 'Schmidt' and 'Seehafer' 
denote the two methods used to compute the initial field for the magnetofrictional extrapolations; 'TD_ex' is the external potential field generated 
by the magnetic sources and dipoles in the TD model. For the unstable case, the reference fields after the initial relaxation (Reference ini) and at 
the end (Reference fin) of the MHD simulation are analyzed. Similarly, two extrapolated fields are reported, corresponding to the maximum (MF 
Schmidt max) and minimum (MF Schmidt min) energy of the reconstructed field fSect. |431 l. All quantities are defined in Sect.|4]and computed in 
the analysis domain. The derivatives involved in ctj and (\fi\) are computed to fourth order accuracy, except for the reference fields and for the 
potential and linear fields, where a second order diff'erencing is used, and one-sided derivatives of corresponding order are employed at, and 
where necessary also just above, the plane {z = 0}. Missing values represent ill-defined quantities for the respective equilibria. 



figurations containing an HFT, a BP, or none of these features. 
The energy, flux rope twist and position, and helicity are ob- 
tained at relative error levels of ^ 1, 3, 4, and ~ 5%, respectively, 
even for a high twist near the threshold of the helical kink insta- 
bility, slightly exceeding one full turn. 

These quantities are found to depend only weakly on the ini- 
tial field chosen, except for the apex height of the flux rope which 
shows a displacement from the true value of up to ^ 20% (if a 
linear force-free field is used). Starting the extrapolation from 
the conventional potential field, which assumes vanishing flux 
outside the magnetogram, yields the most accurate results for 
the considered equilibria. Other potential fields result in some- 
what lower accuracy, and a linear force-free field scores worse 
but still produces a rather acceptable reconstruction. 

The high accuracy is reached using vector magnetograms 
that have a realistic low-flux margin of width ~ 1/4 the size 
of the strong-flux area and are exactly flux balanced. The ef- 
fects that result from violating the latter condition require further 
study. 

Attempting to extrapolate the magnetogram of an unstable 
modified TD equilibrium, the magnetofrictional relaxation pro- 
duces an evolution in two phases of opposing tendency which 
clearly reflect the unstable nature of the equilibrium. A phase of 
rising free energy yielding a partial buildup of the flux rope is 
superseded by a decrease of the energy back to nearly the initial 
(potential-field) value, destroying the flux rope. 
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